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Abstract. This paper develops and analyzes finite element Galerkin and 
spectral Galerkin methods for approximating viscosity solutions of the fully 
nonlinear Monge-Ampere equation det(D 2 u°) = / (> 0) based on the vanish- 
ing moment method which was developed by the authors in 1171 1151 . In this 
approach, the Monge-Ampere equation is approximated by the fourth order 
quasilinear equation — eA 2 u £ + det D 2 u E = f accompanied by appropriate 
boundary conditions. This new approach allows one to construct convergent 
Galerkin numerical methods for the fully nonlinear Monge-Ampere equation 
(and other fully nonlinear second order partial differential equations), a task 
which has been impracticable before. In this paper, we first develop some 
finite clement and spectral Galerkin methods for approximating the solution 
u e of the regularized fourth order problem. We then derive optimal order 
error estimates for the proposed numerical methods. In particular, we track 
explicitly the dependence of the error bounds on the parameter e, for the error 
u £ — u s h . Due to the strong nonlinearity of the underlying equation, the stan- 
dard perturbation argument for error analysis of finite element approximations 
of nonlinear problems does not work here. To overcome the difficulty, we em- 
ploy a fixed point technique which strongly makes use of the stability property 
of the linearized problem and its finite element approximations. Finally, using 
the Aygris finite element method as an example, we present a detailed nu- 
merical study of the rates of convergence in terms of powers of e for the error 
u° — u e h , and numerically examine what is the "best" mesh size h in relation 
to e in order to achieve these rates. 



Fully nonlinear partial differential equations (PDEs) are those equations which 
are nonlinear in the highest order derivative(s) of the unknown function(s). In the 
case of the second order equations, the general form of the fully nonlinear PDEs is 
given by 



where D 2 u°(x) and Du°(x) denote respectively the Hessian and the gradient of u° 
at x e £1 C R". F is assumed to be a nonlinear function in at least one of entries 
of D 2 vP . Fully nonlinear PDEs arise from many scientific and engineering fields 
including differential geometry, optimal control, mass transportation, geostrophic 
fluid, meteorology (cf. [7J |5J [HI HHJ HZ] and the references therein). 
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In this paper we focus our attention on a prototypical fully nonlinear second 
order PDE, the well-known Monge-Ampere equation. Our goal is to develop and 
analyze finite element and spectral Galerkin methods for approximating viscosity 
solutions of the Dirichlet problem for the Monge-Ampere equation (cf . [22 ) : 

(1.2) det(D 2 u°) = f in n C R", 

(1.3) u° = g ondfl, 

where f2 C R™ is a convex domain with smooth or piecewise smooth boundary dfl. 
det(D 2 u°(x)) denotes the determinant of D 2 u°(x). Clearly, the Monge-Ampere 
equation is a special case of (|1.1| ) with F(D 2 u° , Du° ,u° , x) — det(D 2 w°) — /. We 
refer the reader to [4"1 151 12T1 122] and the references therein for the derivation, appli- 
cations, and properties of the Monge-Ampere equation. 

For fully nonlinear second order PDEs, in particular for the Monge-Ampere equa- 
tion, it is well-known that their Dirichlet problems do not have classical solutions 
in general even /, g and dfl are smooth, if f2 is not strictly convex, (see [21]). So it 
is imperative to develop some weak solution theories for these problems. However, 
because of the fully nonlinearity of these equations, unlike in the case of linear and 
quasilinear PDEs, the usual weak solution theory based on integration by parts 
does not work for fully nonlinear PDEs, hence, other "nonstandard" weak solu- 
tion theories must be sought. In the case of the Monge-Ampere equation, the first 
such theory was due to A. D. Aleksandrov, who introduced a notion of generalized 
solutions and proved the Dirichlet problem with / > has a unique generalized 
solution in the class of convex functions (cf. [I], also see [5]). But major advances 



on analysis of problem (1.2 1- ( 1.3 1 has only been achieved many years later after the 
introduction and establishment of the viscosity solution theory (cf. [T2l [22] ) . We 
recall that the notion of viscosity solutions was first introduced by Crandall and Li- 
ons |llj in 1983 for the first order fully nonlinear Hamilton- Jacobi equations. It was 
quickly extended to second order fully nonlinear PDEs, with dramatic consequences 
in the wake of a breakthrough of Jensen's maximum principle [23] and the Ishii's 
discovery [24] that the classical Perron's method could be used to infer existence 
of viscosity solutions. It should be noted that there also exist nonconvex solu- 
tions to problem (|1.2[ )-( 1.3 ), so if the convexity requirement is dropped, solutions 



to problem ( 1.2 1- ( 1.3 > are not unique. We shall refer this uniqueness property as 
conditional uniqueness. Nonuniqucncss or conditional uniqueness is often expected 
for the Dirichlet problems of fully nonlinear second order PDEs. It should also be 
noted that unlike in the case of fully nonlinear first order PDEs, the terminology 
"viscosity solution" loses its original meaning in the case of fully nonlinear second 
order PDEs. 



For the Dirichlet Monge-Ampere problem ( 1.2 )-( 1.3 1 with / > 0, we recall that 
[22] a convex function u° 6 C°(f2) satisfying u° — g on dfl is called a viscosity 
subsolution (resp. viscosity supersolution) of ( |1.2| if for any <p G C 2 there holds 
det(D 2 ip(xo)) < f(xo) (resp. det(D 2 ^(xo)) > /(xo)) provided that u° — ip has a 
local maximum (resp. a local minimum) at xq e 17. «° £ C°(f2) is called a viscosity 
solution if it is both a viscosity subsolution and a viscosity supersolution. Form this 
definition we can see that the notion of viscosity solutions is not variational and 
not defined using the more familiar integration by parts approach which in fact is 
not possible because of the fully nonlinearity of the PDE. On the other hand, the 
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strategy of shifting derivatives onto the test functions at extremal points in the def- 
inition of viscosity solutions can be viewed as a "differentiation by parts" approach. 
It can be shown that a viscosity solution satisfies the PDE in the classical sense at 
every point where the viscosity solution has second order continuous derivatives. 
However, the theory does not tell what equations or relations the viscosity solution 
satisfies at other points. Numerically, this non- variational nature of the viscos- 
ity solution theory poses a daunting challenge for computing viscosity solutions 
because it makes it impossible for one to directly approximate viscosity solutions 
using any Galerkin type numerical methods including finite clement, spectral and 
discontinuous Galerkin methods, which all are based on variational formulations of 
PDEs. In addition, it is extremely difficult (if all possible) to mimic "differentiation 
by parts" approach at the discrete level, so there seems have no hope to develop a 
discrete viscosity solution theory. 

To overcome the above difficulties, recently in |151 117] we introduced a new 
approach, called the vanishing moment method, for establishing the existence of 
viscosity solutions for fully nonlinear second order PDEs, in particular, for prob- 



lem ( 1.2 )- ( 1.3 1 . In addition, this new approach gives arise a new notion of weak 
solutions, called moment solutions, for fully nonlinear second order PDEs. Further- 
more, the vanishing moment method is constructive, so practical and convergent 
numerical methods can be developed based on the approach for computing the vis- 



cosity solutions of fully nonlinear second order PDEs such as problem (1.2l-(1.3l. 
The main idea of the vanishing moment method is to approximate a fully nonlinear 
second order PDE by a quasilinear higher order PDE. The notion of moment solu- 
tions and the vanishing moment method are natural generalizations of the original 
definition of viscosity solutions and the vanishing viscosity method introduced for 
the Hamilton- Jacobi equations in [TT]. We now briefly recall the definitions of mo- 
ment solutions and the vanishing moment method, and refer the reader to |15] 117] 
for a detailed exposition. 

Firstly, the vanishing moment method approximates the fully nonlinear equation 



(111 by the following quasilinear fourth order PDE: 



(1.4) -eA 2 u £ + F(D 2 u e ,Du e ,u e ,x) = (e > 0), 

which holds in domain fl. Suppose the Dirichlet boundary condition u° — g is given 
on dfl, then it is natural to impose the same boundary condition on u £ , hence, 

(1.5) u £ — g on dfl. 



Secondly, in addition to boundary condition (1.5) one more boundary condition 
must be imposed to ensure uniqueness of solutions. In |15] we proposed to use one 
of the following boundary conditions: 

(1.6) Au £ = e, or D 2 u e v -v = e on dn, 

where v stands for the unit outward normal to dfl. Although both boundary 
conditions work well numerically, the first boundary condition is more convenient for 
standard Galerkin type methods such as finite methods, spectral and discontinuous 
Galerkin methods, while the second boundary condition fits mixed finite element 
methods (cf. [El [18]) better. Hence, in this paper we shall use the first boundary 
condition. We also refer the reader to [T7] for the heuristic argument why these 
boundary conditions were chosen in the first place. 
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To sum up, the vanishing moment method consists of approximating second 
order boundary value problem ( |1.3[ )~( |l~Tj ) by fourth order boundary value problem 
( 1.4 1— ( 1.5 ) , (1.6 1. In the case of the Monge-Ampere equation, this then results in 
approximating boundary value problem ( | 1 - 2 [ ) — ( 1.3 1 by the following problem: 



(1.7) 
(1.8) 
(1.9) 



-eAV+det(L>V) = / 



Au £ 



in f2, 
on dQ, 
on dfl. 



It was proved in [TS] that if / > in fl then problem (1.7)-(1.9) has a unique 
solution u £ which is a convex function over f2. Moreover, u £ uniformly converges 
as e — > to the unique viscosity solution of (1.2l-(1.3l. As a byproduct, this 



also shows that (1.2 1— (1.3| has a unique moment solution which coincides with the 
unique viscosity solution. Furthermore, it was proved that there holds the following 
a priori bounds which will be used frequently later in this paper: 

(1.10) |K||^-=0( e -^), \\u s \\ w *,~ =0(-), ||cof(£>V)|U 



o(i) 



for j = 2,3. Where cof(D 2 u £ ) denotes the cofactor matrix of the Hessian D 2 u £ . 
With the help of the vanishing moment method, the difficult task of computing 
the unique convex viscosity solution of the fully nonlinear second order Monge- 



Ampere problem ( 1.2 1— ( 1.3 1, which has multiple solutions (i.e. there are non-convex 
solutions) , is now reduced to a feasible task of computing the unique regular solution 
of the quasilinear fourth order problem ( 1.7 )—( 1.9 ) . In particular, this allows one to 



use and/or adapt the wealthy amount of existing numerical methods, in particular, 



finite clement methods to solve problem (1.2 ) — ( 1.3 1 via problem (1.7 1— (1.91 



The specific goal of this paper is to develop and analyze Galerkin methods for 
approximating the solution of (1.7l-(1.9l in 2-D and 3-D. When deriving error 
estimates of the proposed numerical methods, we are particularly interested in 
obtaining error bounds that show explicit dependence on e. Aygris confirming 
finite element method in 2-D and Legendre spectral Galerkin method in both 2-D 
and 3-D are specifically considered in the paper although our analysis applies to any 
conforming Galerkin method for problem (1.7l-(1.9l. We note that finite element 



approximations of fourth order PDEs, in particular, the biharmonic equation, were 
carried out extensively in 1970's in the two-dimensional case (see [TO] and the 
references therein), and have attracted renewed interests lately for generalizing the 
well-know 2-D finite elements to the 3-D case (cf. [34] [35] [33] ) and for developing 
discontinuous Galerkin methods in all dimensions (cf. 16, 28). Clearly, all these 



methods can be readily adapted to discretize problem (1.7l-(1.9l although their 



convergence analysis do not come easy due to the strong nonlinearity of the PDE 
(1.7 1. Also, the standard perturbation technique for deriving error estimates for 



numerical approximations of mildly nonlinear problems does not work for problem 
(1.7 1— ( 1.9 ) . We refer the reader to [18, 29J for further discussions in this direction. 



We also like to mention that a few results on numerical approximations of the 
Monge-Ampere equation as well as related equations have recently been reported 
in the literature. Oliker and Prussner [31] constructed a finite difference scheme for 
computing Aleksandrov measure induced by D 2 u in 2-D and obtained the solution u 



of problem ( 1.7 H 1.9 1 as a by-product. Baginski and Whitaker [2J proposed a finite 
difference scheme for Gauss curvature equation (cf. [T7] and the references therein) 
in 2-D by mimicking the unique continuation method (used to prove existence of 



GALERKIN METHODS FOR MONGE-AMPERE EQUATIONS 



5 



the PDE) at the discrete level. In a series of papers (cf. [13] and the references 
therein) Dean and Glowinski proposed an augmented Lagrange multiplier method 



and a least squares method for problem (1.7 1— ( 1.9 ) and the Pucci's equation (cf. 



[71 [ST]) in 2-D by treating the Monge- Ampere equation and Pucci's equation as 
a constraint and using a variational criterion to select a particular solution. Very 
recently, Oberman [30 constructed some wide stencil finite difference scheme which 
fulfill the convergence criterion established by Barles and Souganidis in [3] for finite 
difference approximations of fully nonlinear second order PDEs. Consequently, 
the convergence of the proposed wide stencil finite difference scheme immediately 
follows from the general convergence framework of |3J. Numerical experiments 
results were reported in [3TJ |30j El [13], however, convergence analysis was not 
addressed except in [50] . 

The remainder of this paper is organized as follows. In Section [2] we first derive 



the weak formulation for problem ( 1.7 1-( 1.9 1 and then present our confirming finite 



clement and spectral Galerkin methods based on this weak formulation. In Section 
[3] we study the linearization of problem (1.7 1 — ( 1 .0 ) and its Galerkin approxima- 



tions. The results of this section, which are of independent interests in themselves, 
will play an important role in our error analysis for the numerical method intro- 
duced in Section [2] In Section [4j we establish optimal order error estimates in 
the energy norm for the proposed confirming finite element and spectral Galerkin 
methods. Our main ideas are to use a fixed point technique and to make strong use 
of the stability property of the linearized problem which is analyzed in Section [3] 
In addition, we derive the optimal order error estimate in the L -norm for u e — u £ h 
using a duality argument. Finally, in Section [5] we first run some numerical tests 
to validate our theoretical error estimate results. We then present a detailed com- 
putational study for determining the "best" choice of mesh size h in terms of e in 
order to achieve the optimal rates of convergence and for estimating the rates of 
convergence for both u° — u s h and u° — u £ in terms of powers of e. 

2. Formulation of Galerkin Methods 

Standard space notations are adopted in this paper, we refer to [2TJ [TU] for 
their exact definitions. In addition, f2 denotes a bounded domain in R n . (•, •) and 
(■,•) denote the L 2 -inner products on £1 and on dfl, respectively. C is used to 
denote a generic e-independent positive constant. We also introduce the following 
special space notation: 

V := H 2 (fl), V Q := H 2 (fl) R HfcSl), V g := {v € V; v\ m = g}. 



(2.1) 



Testing (1.7) with v € Vq yields 

-s f Au £ Avdx + [ det(D 2 u 6 )vdx = [fvdx - [ e 2 ^"-ds. 
Jn Jn Jn Jon ®v 



Based on (2.1 1, we define the variational formulation of (1.7 1 — ( 1.9 ) as follows: Find 



u e £ V g such that 

(2.2) -e(Au E ,Av) + (det(D 2 u E ),v) = (f,v)-(e 2 ,^) Vv£V . 

\ ° v I dn 

a 

Remark 2.1. We note that det(D 2 ir) = ^<$> e D 2 u £ = A ^ ^ £ ijU XzX] j = 1, 2, ...n., 

i=l 

where <I> £ is the cofactor matrix of D 2 u e . Thus, using the divergence free property 
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of the co] 'actor matrix <I> £ (cf. Lemma 3.1) we can define the following alternative 
variational formulation for (12. 21) . 



(2.3) -e(Au £ ,A^--(^Du s ,D^) = {f,^)-(s 2 ,^-) ^eV . 
n \ av I an 

However, we shall not use the above weak formulation in this paper although it is 

interesting to compare Galerkin methods based on the two different but equivalent 

weak formulations. 



In this paper, we shall consider two types of Galerkin approximations for (2.2 1. 
The first type is the confirming finite element Galerkin method in 2-D. Aygris fi- 
nite element will be used as a specific example of this class of methods although 
our analysis is applicable to other confirming finite element such as Bell clement, 
Bogner-Fox-Schmit element, and Hsich Clough Tocher clement (cf. [10J). The 
second type of the methods is the spectral Galerkin method in 2-D and 3-D. Lc- 
gendre spectral Galerkin method will be used as a specific example although our 
analysis also applies to other spectral methods (with necessary assumptions on the 
domain in the case of Fourier spectral method) . 

To formulate the finite clement method in 2-D, let Th be a quasiuniform trian- 
gular or rectangular partition of Q, C R 2 with mesh size h £ (0, 1) and V h C V 
denote a confirming finite element space consisting of piecewise polynomial func- 
tions of degree r(> 5) such that for any v £ V D H s (£l) 

(2.4) inf \\v-v h \\ Hi <^- j '|H|h», j = 0,l,2;^ = min{r+l,s}. 

v h <£V h 

We recall that r — 5 in the case of Aygris element (cf. [TO])- 
Let 

Vl l = {v h e V h ; v h \ dn = g}, V h = {v h £ V h ; v h \ an = 0}. 

Based on the weak formulation ( |2.2| , we define our finite element Galerkin method 
as follows. Find ul £ Vj} such that 



Vv h e V h . 



(2.5) - e(Aui,Av h ) + {det(D 2 u%),v h ) = (f,v h ) ~ (e 2 ,^) 

\ ov I an 

To formulate the spectral Galerkin method, we assume that £1 is a rectangular 
domain and let Lj denote the jth order Legendre polynomial of single variable and 
define the following finite dimensional spaces: For iV > 2, let 

V := span{Lo(^i), ^2(^1), • • ■ , ^jv(^i)} when n = 1, 

V N := span{Li(xi)Lj(x2); i, j = 1, 2, • • • , N} when n = 2, 

V N := sp&n{Li(xi)Lj(x2)L).(x3); i, j, k = 1, 2, • • • , N} when n = 3. 

It is well-known that V N has the following approximation property (cf. [S]): 

(2.6) inf \\v- v N \\ H 3 < CW-'WvWhs, < j < minis, N + 1} 

v N <£V N 

for any v £ V H s (£l). 



Again, based on the weak formulation (2.2), our spectral Galerkin method for 



(|2.2| is defined as seeking u £ N £ V N n V g such that for any v N £ V N n V 



(2.7) - e(A^, Av N ) + (det(D 2 u%),v N ) = (/, v N ) - (e 2 , 



<><> 
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Clearly, Galerkin methods (2.51 and (2.7 1 have the exact same form as the vari- 
ational problem (2.2 1. The only difference is that the infinite dimensional space V 
in (2.2 1 is replaced by the finite dimensional subspace V h and V N , respectively. 
Letting h := -4, u h := u N and V h := V N in the definition of the spectral method, 
(2.7 1 can be rewritten as the exactly same form as (2.5 1. For this reason, we shall 
abuse the notation by using u h to denote the solution of (2.5 1 and the solution of 
(2.7 1 with understanding that h = j^. Since the convergence analyses for (2.51 and 
for (2.7) are essentially same, we shall only present the detailed analysis for (2.5) 
and make comments about that for (2.7 1 in case there is a meaningful difference. 

Let u e be the solution to (2.2) and u\ the solution of (2.51 or (2.7 1. As men- 
tioned in Section [l] the main task of this paper is to derive optimal error estimates 
for u e — it?. To this end, we first need to prove existence and uniqueness of u e h . 
It turns out both tasks are not easy due to the strong nonlinearity in (2.5 1 and 
( |2.7| . Unlike the continuous PDE case where u e is proved to be convex for all e (cf . 
|17j ). it is not clear whether u\ preserves the convexity even for small h. Without 
a guarantee of convexity for w|, we could not establish any stability result for u e h . 
This is the main obstacle for proving existence and uniqueness for (2.5) and (2.7 1. 
In addition, again due to the strong nonlinearity, the standard perturbation tech- 
nique for deriving error estimate for numerical approximations of mildly nonlinear 
problems does not work here. To overcome the difficulty, our idea is to adopt a 
combined fixed point and linearization technique which was used by the authors in 
|19j . where a nonlinear singular second order problem known as the inverse mean 
curvature flow was studied. We note that by using this technique we are able to 
simultaneously prove existence and uniqueness for u E h and also derive the desired 
error estimates. In the next two sections, we shall give the detailed account of this 
technique and apply it to problems (2.5 1 and (2.7 1. 



3. Linearization and its Finite Element Approximation 



To analyze (2.5l and (2.7l, we shall study the linearization of (1.7) to establish 
the required technical tools. First, we recall the following divergence-free row prop- 
erty for the cofactor matrices, which will be used frequently in later sections. We 
refer the reader to [HJ p. 440] for a short proof of the lemma. 

Lemma 3.1. Given a vector-valued function v = [v\,V2, ■ ■ ■ ,v n ) : f2 — > M. n . As- 
sume v € [C 2 (f2)] n . Then the cofactor matrix cof(Dv) of the gradient matrix Dv 
o/v satisfies the following row divergence-free property: 

n 

(3.1) div(cof(£>v))i = 9 Xj (cat(pv))ij = for i = 1, 2, • • • , n, 

where (cof(-Dv)), and (cof(Dv))jj denote respectively the ith row and the (i, j)-entry 
o/cof(Dv). 

3.1. Linearization. It is easy to check that for a given smooth function w there 
holds 

(3.2) det(L>V + tw) = det(L>V) + ttr(t> £ D 2 w) + ■■■ + t n det(D 2 w). 
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The linearization of M £ (u £ ) := — eA 2 u e + det(D 2 u e ) at the solution u £ is given by 

M £ (u £ + tw) -M e (u e ) 



L u e(w) : 



lim 

f->0 



-eA 2 w + § £ : D 2 i 



-eA 2 w + div($ £ Dw), 



where <I> e denotes the cofactor matrix of D u £ and we have used Lemma 3.1 to get 
the last equality. Also, using Lemma 3.1 it is easy to check that L u e is self-adjoint, 
i.e., the adjoint operator L* E of L u e coincides with L u e. 
We now consider the following linear problem. 

(3.3) L u e(v) — tp in ft, 

(3.4) v = on dfl, 

(3.5) Av = q on dfl. 

Multiplying the PDE by a test function w € Vq and integrating over we get the 



following weak formulation for ( 3.3 1— ( 3.5 1 : Find v £ Vq such that 

/ dw \ 

(3.6) B[v,w} = ( V ,w)-e(q, — ) VweV 0} 

where 



B[v,w] ■= £ AvAwdx+ / $ £ DvDi 
Jn io 



>w dx. 

in Jn 

The next theorem ensures the well-posedness of the above variational problem. 



Theorem 3.2. Suppose d£l € C 2 . Then for every ip 6 V * and q € H~^{d£l) 
there exists a unique solution v Vq to problem (3.6 1. Furthermore, there exists a 
positive constant Ci(e) such that 



(3.7) 



Miff 2 



<Ci(s)(\ 



Proof. It is easy to check that -B[-,-] is a bounded and coercive bilinear form on 
Vq x Vq with coercive constant C^e) = 0(e). Also, the right-hand side of (3.6) 



defines a bounded linear functional on Vq (see [29] for details). So the assertions 
of the theorem follows immediately from an application of Lax-Milgram Theorem 

HB. □ 

For more regular data, the above theorem can be improved to the following one 
(see for a similar proof). 

Theorem 3.3. Suppose ip £ (H s (fl) n Hq(CI))*, q e H s ~i (dty, dtt £ C s , (s > 2) 
and v is the unique solution to (3.6 I. Then v £ H s (fl) n Hq(Q,), and there exists a 
C s (e) > such that 



(3.8) 



< c s{z) {\W\\{H°nHl(n)y + \\q\\ h . 



■a (an) 



Remark 3.1. From the definition ofC2(s), we see that in Theorem 3.2 Ci(e) = 

remains 



3.3 



O(i). Currently, the explicit dependence of C s (e) on e in Theorem 
unknown for generals > 3. Finally, we note that Theorem \3.S\ can easily be extended 
to the case of nonhomogeneous boundary data which we summarize in the following 
theorem. 
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Theorem 3.4. For every tp 6 V* , g £ H^^dfl), q £ H ^(dfl), there exists a 
unique weak solution w £ V to 

L u e(w) — p in f2, 

w = g on <9f2, 

Aw = q on d£l. 
Furthermore, there exists C > swc/i that 
C 



(3.9) 



M^ 2 < 



M(Hi l nmy + \\g\\ H i u 



(on) 



We refer the reader to [23j for a similar proof. 



3.2. Finite element approximation of linearized problem. Let V h C Vq be 

one of the finite-dimensional subspace of Vq as defined in the previous subsection 
(e.g. Aygris finite element and Legendre spectral element), and v £ Vq denote the 
solution of (3.6 1. Based on the variational equation (3.6 1, our Galerkin method for 
(3.3 1 is defined as seeking Vh £ V such that 

dw h 



(3.10) 



B[v h , w h ] = (cp, w h ) - £ ( q, 



on 



Our objective in this subsection is to first prove existence and uniqueness for 
problem (3.101 and then derive optimal order error estimates in various norms. 

Theorem 3.5. Suppose v £ VoC)H s (il) (s > 3). Then there exists a unique Vh £ Vq 1 
satisfying (3.101. Furthermore, we have the following estimates: 

(3-11) Kll*'(n) < C 3 (e) (|| 

<P\\ (H^nH 2 )* + \\q\\ H -han) ) ' 

(3.12) ||v - v h \\ HHQ ) < C 4 (e)h e - 2 \\v\\ Ht{n) , 

(3.13) ||v - v h \\ H t in) < Cs{e)h l - x \\v\\ H z {n) , 

(3.14) \\v-v h \\ L 2 W < C 6 (e)h e \\v\\ 

where i = min{r + l,s} in the case of the finite element Galerkin method, and 
I = min{./V + 1, s} in the case of the spectral Galerkin method. 



Proof. Estimate (3.11) follows immediately from setting Wh = in (3.10) and 
using the coercivity of the bilinear form B[-, ■]. 

To derive the error estimate in the i? 2 -norm, we use the error equation, 

B[v-v h ,w h } = Vw h eV h . 

Using the coercivity property of B[-, •], we have 

(3.15) C 2 (e)||« - v h \\ 2 H 2 < B[v -v h ,v- v h ] = B[v - v h ,v] - B[v - v h ,v h ] 

= B[v - v h ,v] = B[v -Vh,v- Ihv], 

where C 2 (e) = 0(e), IhV denotes the Galerkin interpolant of v onto Vq. On noting 
that 

C 

B(v -v h ,v- I h v) < -^\\v - v h \\ H 2\\v - hv\\ H 2, 



we have 



v - V h \\ H 2 



< 



C 



\\v-I h v\\ H 2 <C 4 (e)/i z - 2 ||z;|| ffi , 
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where C 4 (e) = O [e~*\ . Thus, ( [3.12| holds. 

Next, we derive the i/^-norm error estmate using a duality argument. Define 
e/s := v — Vh and consider the following problem: 





Ae ft 


in f2, 


il> = 





on 


Atp = 





on d£l. 



Using (3.8 1, we have 

\\i>\\ H »<C a (e)\\Ae h \\ H -x. 
Since ||Ae^|| ff -i = sup{(Ae/ l; u)\ u S Ho(fi), < 1}, we have 

(Ae h ,u) = (Ve hs V«) < || V e/l || L 2 1| Vw|| L 2 < ||Ve h .|| L 2|| u || H i = ||Ve & ||joa. 

It follows that 

||Ae h || H -i < HVefcHia. 

Thus, 

||Ve ft ||| 2 = (Ae h ,e h ) = B[e h ,ip] = B[e h , tp - I h ip] 
C hC 

< /l %M|j Ae;i || ff _ 1 || e , i || ff2 < hC ^\\Ve h \\ L2 \\e h \\ H2 . 



Hence, 



\\Ve h \\ L * < /lC ^||e /t || ff2 . 



Combining the above with ( |3.12[ ), we get |3T3| ) with C 5 = O (C s (e)£- 2 ) . 
To derive the error in the L 2 -norm, we consider the following problem: 





= e h 


in f2, 




= 


on dfl, 


Aip 


= 


on dfl. 



On noting (3.8 1 implies that 



H\\ H * < C s (e)\\e h \ 



L 2 ■ 



Thus, 



\\ e h\\ 2 L 2 = {e h ,V- v h ) = B[v - Vh,ip] = B[v -v h ,ip- I h ip] 



C 

< —\\v-V h \\ H 2\\%p-I h lp\\ H 2 

y/e 

h 2 C u . h 2 CJe) ,., 

<-r\\v-v h \\ H 2\\il>\\ H *< sV '- 



Dividing by || e/j. || ^2 , we get (3.141 with Cq = O (C s (e)e 2 ). The proof is complete. 



□ 



Remark (a) In the case of Aygris finite element method, r = 5 in (3.11 |-(3.14l. 

(b) In the case of Legendre spectral method, h = -h in (3.11 1-(3.14|, where N 
stands for the highest degree of polynomials in V N . 
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4. Error Analysis for Galerkin Methods (2.5) and (2.7 1 



The goal of this section is to derive optimal order error estimates for the Galerkin 
methods (2.51 and (2.7). Our main idea is to use a combined fixed point and 
linearization technique which was introduced by the authors in |19j . Once again, we 
only present the detailed analysis for (2.51 since the analysis for (2.7 1 is essentially 
same. 

First, we define a linear operator T : V g — > V g . For any S V g , let T(w/ l ) S 
V g h denote the solution of following problem: 

(4.1) B[w h - T(w h ),fa] = e(Aw h ,Afa) - (det(D 2 w h ), fa) 



dfa 
dv 



Vfa e v n 



mi 



By Theorem 3.5 we see T(wh) is uniquely defined. Notice that the right-hand 
side of (4.1 ) is the residual of to equation (2.5 1. It is easy to see that any fixed 
point Wh of the mapping T (i.e. T(wh) = Wh) is a solution to problem (2.5 1 and 
vice-versa. In the following we shall show that the mapping T indeed has a unique 
fixed point in a small neighborhood of IhU £ . To this end, we set 

B h {p) := {v h G V g h ; \\v h - I h u e \\ H 2 < p}. 

In the rest of the section, we assume u £ G H s (£l) and set I = min{r + 1, s}. 

Lemma 4.1. There exists a constant CV(e) = O (e~ n ) > (n = 2, 3) such that 

(4.2) \\I h u e - T{I h u s )\\ m < C 7 {e)h l ' 2 \\u e \\ Hl . 

Proof. To simplify notation, let lu^ := I^u 6 — T(// l u 6 ) and let rj^ :— I^u 6 — u £ . We 
then have 

duih 



B[uj h ,uj h ] = e(A(I h u e ),Auj h ) - (det(D 2 (I h u E )) -f,uj h )-(e 



ch 



on 



= e(Arj h , Au> h ) + (det(£>V) - det{D 2 (I h u e )) , u) h ) 
< ellATjfcHiallAwfcll^ + \\det(D 2 u £ ) - det(D 2 (I h u e ))\\ L 2\\u;\\ L 2 
Using the Mean Value Theorem we have 

det(D 2 (I h u e )) - det(Z?V) = $ : D 2 r, hl 

where <l = cof(TD 2 (I h u E ) + (1 - r)D 2 u e ) for some r G [0, 1]. 
Next, when n = 2, we bound as follows: 

||<I £ || L ~ = ||cof(r D 2 (I h u e ) + (1 - T)D 2 u e )\\ L ~ 
= \\TD 2 (I h u £ ) + (l-T)D 2 u e \\ L ~ 
< \\D 2 (I h u')\\ L ^ + \\D 2 u'\\ L ^ 
C 



< C\\D 2 u s \\ 



< 



Similarly, when n = 3, we can show that = 0(e~ 2 )- Hence, 

C 

B[uj h ,w h ] < s\\Ar) h \\ L 2\\Au; h \\ L 2 + ^^jll-D^IU 2 IkftlU 2 
C 

- ^tII^II^II^IIh 2 - 
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Using the coercivity of the bilinear form £?[•,•] we get 



Thus, (|42j holds with (7 7 (e) = e „_, c Ci(E) = O (£~ n ). □ 

Lemma 4.2. There exists an ho > and < po < 1 suc/i £/ia£ /or h < /io, £/ie 
mapping T is a contracting mapping in the ball Bh(po) with a contraction factor g. 
That is, for any € Bh(po), there holds 

(4.3) ||T(« h ) - TK)|| ff2 < i||« h - ^|| ff2 . 

Proof. For any -0^ G Vq , using the definition of T(vh) and T(wh) we get 

B[T(« h ) - T(wh),iph] = {® £ D{v h - w h ),D^ h ) + (det(D 2 v h ) - det(D 2 w h ),^ h ) . 

Let rf, iv? denote the standard mollifications of Vh and Wh, respectively. Adding 
and subtracting these terms and using the Mean Value Theorem yield 

B[T(v h ) - T(wh),i>h] 

= ($ £ (Dv h - Dw h ),D^ h ) + {dct(D 2 v h ) - det(D 2 w h ),tp h ) 
= ($ s (Dv h - Dw h ),D^ h ) + (det(£>X) - det(£> 2 <).^) 

+ (det(D 2 v h ) - det(L> 2 <),^) + (det(£> 2 <) - det(D 2 w h ), $ h ) 
= (t> e (Dv h - Dw h ),D^j h ) + : " ^ 2 <): WO 

+ (det(£>V) - det(L> 2 0, V/,) + (det(£>X) - det(D 2 w h ), i() h ), 

where # h = cof(D 2 u£ + t(D 2 w£ - D 2 v^)), r € [0, 1]. 
Using Lemma|3.1| we have 



B[T(v h )-T(w h ),i> h ] 

= - V h )(Dv h - Dw h ),Dt/j h ) + (V h (Dv h - Dv%), Diph) 
+ (* h {Dw£ - Dw h ),ip h ) + (det(D 2 v h ) - det(D 2 v^)^ h ) 
+ (det(L>X) - det(D 2 w h ), i/> h ) 

< C*{||$ £ - * fc ||ia||«h - W h \\ H 2\\lp h \\ H 2 + W^hW^UhWm [IK - v£\\ B * 

+ IK - KWhA + [||det(L> 2 ^) - det( J D 2 <)|| L 2 

+ ||det(£>V) - det(U 2 0|| z >] ||Vh|U»}, 

where we have used Sobolev's inequality. 

Next, we derive an upper bound for ||<I> e — ^Uz^ when n=2 as follows: 

W - * fc || ia = \\& - cof(D 2 < + t(D 2 w£ - D 2 v£))\\ L 2 

= ||DV - (D%£ + r(U 2 < - ^ 2 <))|| i2 

< ||£>V - D 2 (I h u e )\\ L 2 + \\D 2 (I h u £ ) - L>V|| i2 

+ 2\\D 2 v h - D 2 v£\\ L 2 + \\D 2 w h - D 2 w£\\ L 2 + \\D 2 v h - D 2 w h \ 

< Ch e - 2 \\u e \\ H i + 2 Po + 2\\v h - v%\\ H 2 + |K - <||j?2. 
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Similarly, when n=3, we can show 

- *h\\ L * < j(h e - 2 \\u e \\ H t + 2 Po + 2\\v h - v£\\ H2 + \\w h - <|| ff2 ) 

Using this result above, we have 

B[T(v h )-T{w h ),iP h ] 
C 



< 



£ n-2 



2 \\ u£ \\h* + Po + IK - KWm + \\wh - WhWh 2 ] IK - W h \\ H 2 
+ \\*h\\u> [\\v h - v£\\m + IK - <\\h*] 

+ [||dct(^ 2 %) -dctp 2 <)llL 2 + lldetp 2 ^) -detpX)IU 2 ]}llV^II^- 
Setting p — > yields 

B[T{v h -T{w h ),ip h ] <^(h e - 2 \\u e \\ H e+ Po ) \\v h -w h \\ m \\i> h \\ H 2. 
Using the coercivity of the bilinear form £?[■, •] we get 



C 



\\T(v h ) - T(w h )\\ H 2 < {h e - 2 \\u^\\ He + p ) \\v h - w h \\ 



H 2 ■ 



C 2 (e)e«- 2 

Choosing p — c ' 2 ^ and ho = ( 7jpl| e ) ' 2 ■, then for h < h Q we have 

\\T(v h ) ~ T(w h )\\ H 2 < ~\\v h - w h \\ H 2. 
The proof is complete. □ 
Theorem 4.3. Let pi — 2CT{e)h l ^ 2 \\u e \\ui . Then there exists an hi > such that 



for h < min{/i :^i}; there exists a unique solution uf of (2.5) in the ball Bh(pi). 
Moreover, there exists a constant Cs(e) = O (s~ n ) such that 

(4.4) \\u £ -ut\\ H 2 <C 8 (e)h e - 2 \\u £ \\ H z, t = min{r + 1, s}. 

2n — l 1 

Proof. Let Vh <E Bh(pi) and hi = ( tc^m ~ ) e ~ 2 ■ Then h < min{/io, fei} implies that 
Pi 5: Po- Thus, using the triangle inequality and Lemmas |4.1| and |4.2| we have 



\\I h u' - T{v h )\\ H 2 < \\I h u" - T(I h u")\\ H 2 + \\T(I h u s ) - T(v h )\\ H 2 

< C 7 (e)h e - 2 \\u\\ Hf . + hhu^ -v h \\ H , < £ + § = Pl . 



Hence, T(vh) S Bh{pi). In addition, from Lemma 4.2 we know that T is a con 



tracting mapping. Thus, the Brouwer fixed Theorem [21 infers that T has a unique 



fixed point u\ G Bh(pi), which is the unique solution to (2.5) 



To get the error estimate, we use the triangle inequality to get 

\\u e - ul\\ H 2 < \\u e - I h u e \\ H 2 + \\I h u e - u £ h \\ 

< Ch l - 2 \\u\\ Hl + Pl = C s (e)h*- 2 \\u\\ Hi , 

where C$(e) := CCV(e) = O The proof is complete. □ 
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Theorem 4.4. In addition to the hypothesis of Theorem \4-3[ assume that the 
linearized equation is H -regular with the regularity constant C s {e). Then there 
holds 

(4.5) \\u £ -ut\\ L * < C 9 (e)A\u s \\ Ht + e 2 - n C 8 (e)h 2 ^\\u £ \\ 2 He 

where Cg(e) = C s (e)Cs(e). 

Proof. Let e| := u e — u e h and it? denote a standard mollification of u e h . It is easy 
to verify that e £ h satisfies the following error equation: 

£ (Ae|,A^) + (detpVj -detpV),^) = VV> ft € ttf. 



Using the Mean Value Theorem and Lemma |3.1| we have 

- e(Ae%,A^ h ) + (det(D 2 u £ f) - det(D 2 u £ ), i/> h ) + (det(D 2 u £ h ) - det(D 2 u £ ^ l ),^ h ) 
= e(Ae|, Ai/> h ) - ($D(u £ h <* - u £ ),D^ h ) + (det(D 2 u £ h ) - det(D 2 u £ ^), ij> h ), 

where $ = cof(£> 2 <' M + r(L>V - D 2 u^)), r € [0, 1]. 

Next, The H 4 -regular assumption implies that (cf. Theorem 3.2 1 there exists a 
unique solution to the following problem: 



A.« WO 


= e| 


in f2, 




= 


on <9f2, 


A?/> 


= 


on d£l. 


Ml** 


< C s {e)\ 





Moreover, there holds 
(4.6) 
Thus, 

\K\\h = = e(AelA^) + (^ £ Dip, De\) 

= s(Ae e h , A{i> - I h tf,)) + ($ e De%, D(tf> - I h r/>) + e{Ae%, A{I h ^)) 

+ ($ e De e h ,D(I h ip)) - e{Ae%, A{I h ^ h )) - {$>D{u £ - u £ ^),D(I h ^)) 

- (detpVj - det{D 2 u £ ^),I h ^) 
= e(Ae%, A{i> - I h ^)) + ($ £ De s h , Dty - I h ^)) 

+ (<S> £ De £ h - $D(u £ - u£"), -0(4 V)) " (det(D 2 u £ h ) - det(£> 2 <'"), I h ^) 
= s(Ae £ h , A(V> - h^)) + ($ e Del, - I h ^)) + ((<T - $>)De £ h , D{I h $)) 

+ ($-D«' M - u £ h ),D(I h ^)) + (det(I?V/) - det(DX). W 

< ||Ae||U 2 ||AW-4V)IU= + C||$ s |U2|| e ||| H2 ||^-/^|U 2 

+ C||0> e - $|| i2 ||e^|| H2 ||J^IU= + - ui\\ H 4hTp\\H* 

+ ||det(D a O -det^X'^MI-MI^, 
where we have used Sobolev's inequality. Continuing, we have 

(4.7) ||eJJ| 2 2 < c{sh 2 \\e £ h \\ H i + h 2 \\$ £ \\ L . \\e £ h \\ H , + ||<F - *|U»|KIU> 

+ ||*|| L 2||<-" - u £ h \\ L 2 + ||det(D 2 <) - dBtpVjHiaJllVHfl*. 
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We now bound ||<I> £ — $||l 2 for n—2 as follows: 
(4.8) ||$ £ - $|| L 2 = ||cof(£>V) - cof(D 2 u e ^ + t(L>V - D 2 u £ ^ l ))\\ L 2 



)V - D 2 u e ^ + T{D 2 u\f - £>V) 

l2„,e r>2„,e| _ , || n2„,e n 2„ £ .M| 



< ||DV - £>Xlb + IPX - ^X'lb 

+ \\D 2 u^ - D 2 ui\\ L2 + \\D 2 u% - DV|| ia 
= 2||£>V - ^ 2 <|| i2 + 2||£> 2 < - ^H'HIl- 
Similarly, when n=3, we have 



(4.9) 



Substituting (4.8 )-(4.9 1 into (4.7 1 we obtain 



ley 2 ,, < c(eh 2 \\el\\ H , + h 2 \m\ L * \\e%\\ Ha + ^{hhW* + \\u% u^\\ H .)\\el\\ H 2 
+ m^Wu^ -ul\\ H 2 + ||det(D 2 <) - det(I? 2 «^)|U 3 ) 



Setting /i — > and using (4.6) yield 

\K\\h < C (eh 2 \\e%\\ H , + h^ L qei\\ H *+e 2 - n \\ei\\ 2 H2 ) 

< C s (e) (eh 2 \\el\\ H2 + h 2 \m\ L ,\\el\\ H , + s 2 - n \\e%f H2 ) \\el\\ L ,. 



Hence, 



V < C s (e)(eh 2 \\el\\ H2 + h 2 \\*^„\\e%\\ H 2 + e 2 - n \\e E h \\ 2 Ha ) 



-2-ri|| e ||2 
\ e h\\H 2 



e h\\Hi + - 

\\ut\\ He +e 2 - n C 8 (e) 2 h 2 ^\\uC\\ 2 He 



C 9 (e) — \\u e \\ Hl + e 2 - n C 8 (e)h 2 ^\\ u c\\ Ht 



where C 9 := C s (e)C 8 (e). 



a 



Remark It can be shown that the corresponding error estimates for spectral 
Galerkin method (2.7 1 are 

(4.10) |K - u s N \\ H2 < C 8 {e)N 2 - l \\v?\\ Ht , 



(4.11) 



\\u e -ul\\ L 2 <C 9 (e) iV-V5||^||^+ £ 2 -"C 8 ( £ )7V 4 - 2£ ||^|| 2 



H< 



provided that u € £ H s (£l). Where £ = minjiV + 1, s}. We refer the reader to [25] 
for a detailed proof. 

5. Numerical Experiments and Rates of Convergence 

In this section, we provide several 2-D numerical experiments to gauge the ef- 
ficiency of the finite element method developed in the previous sections. We nu- 
merically find the "best" choice of the mesh size h in terms of e, and rates of 
convergence for both u° — u e and u e — u e h . All tests given below are done on the 
domain Q, = [0, l] 2 . We refer the reader to [T7, 29 for more extensive 2-D and 3-D 
numerical simulations. 
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£ 


lluf - M U || r! 
II h ^ II L z 


llltf - li U ||„-l 
II u 'ft II H 1 


II "'ft Lh II H z 


0.75 


0.109045862 


0.528560309 


3.39800721 


0.5 


0.113340196 


0.548262711 


3.524120741 


0.1 


0.08043631 


0.401646611 


3.071852861 


0.075 


0.06932532 


0.352221521 


2.900677852 


0.05 


0.053925875 


0.283684684 


2.657326288 


0.025 


0.032202484 


0.18559903 


2.270039867 


0.0125 


0.017972835 


0.117524466 


1.928506935 


0.005 


0.007871272 


0.062721607 


1.544066061 


0.0025 


0.004115832 


0.038522721 


1.301171395 


0.00125 


0.002124611 


0.023464656 


1.095145652 


0.0005 


0.00087474 


0.012073603 


0.871227869 



Table 1. Test la: Change of \\vP — u E h \\ w.r.t. e (h — 0.009) 



e 


\Wk-""\\ L 2 


IW-« u ll H i 




£ 


s/i 




0.75 


0.145394483 


0.610328873 


3.651396376 


0.5 


0.226680392 


0.775360561 


4.19090946 


0.1 


0.804363102 


1.270118105 


5.462612693 


0.075 


0.924337601 


1.286131149 


5.542863492 


0.05 


1.078517501 


1.268676476 


5.619560909 


0.025 


1.288099375 


1.173831334 


5.708848032 


0.0125 


1.437826805 


1.051170776 


5.767580991 


0.005 


1.574254363 


0.887017474 


5.806619604 


0.0025 


1.646332652 


0.770454411 


5.819015381 


0.00125 


1.699688442 


0.663680706 


5.82430863 


0.0005 


1.749480266 


0.53994796 


5.826251909 



TABLE 2. Test la: Change of ||u° - u||| w.r.t. e (h = 0.009) 



Test 1: For this test, we calculate — itf || for fixed h — 0.009, while varying e in 
order to approximate ||u e — We use Argyris elements and set to solve problem 



(2.51 with the following test functions 

(a) . u° = e^+f 2 )/ 2 , / = (1 + x 2 + y 2 )e^+^/ 2 , g = e ^ + ^' 2 . 

(b) . u° = x A + y 2 , f = 24a; 2 , g = x* + y 2 . 

After having computed the error, we divide by various powers of e to estimate the 
rate at which each norm converges. Tables [2] and |4] clearly show that ||tt° — u^||jj2 = 
0(ei). Since we have fixed h very small, then ||w° — u £ \\h2 w ||tt° — uf l \\n2 = 0(ei). 
Based on this heuristic argument, we predict that ||u° — it £ ||^2 = O(ei). Similarly, 
from Tables [2] and |i] we see that - u e \\ L 2 w O(e) and - u £ \\ H i w 0(ei). 

Test 2. The purpose of this test is to calculate the rate of convergence of | \u e — u e h \ \ 
for fixed e in various norms. As in Test 1, we use Argyris elements and solve problem 



(2.51 with boundary condition Au E = e on 951 being replaced by Au E = <jf on dfl. 



We use the following test functions: 
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£ 


II \\L Z 


\\ul _ u"\\„l 


lllif — II jjl 
II '"h \\H Z 


0.75 


0.179911089 


0.896016741 


5.98759668 


0.5 


0.177287901 


0.883816723 


5.982088348 


0.1 


0.102586549 


0.549713562 


4.822739159 


0.075 


0.085457592 


0.47264786 


4.537189438 


0.05 


0.063960926 


0.374513017 


4.150185418 


0.025 


0.036755952 


0.24464464 


3.552006757 


0.0125 


0.020291198 


0.157714933 


3.032842066 


0.005 


0.008967657 


0.087384209 


2.451390014 


0.0025 


0.004761813 


0.055425626 


2.080704688 


0.00125 


0.002502224 


0.034885527 


1.762183589 


0.0005 


0.001054596 


0.018689724 


1.410593138 


0.00025 


0.000544002 


0.011565172 


1.189359491 


0.000125 


0.000279021 


0.007112 


0.999863491 


0.00005 


0.000114659 


0.003700268 


0.787092117 



TABLE 3. Test lb: Change of ||u° - u%\\ w.r.t. e (h = 0.009) 



e 


\W%-u u \\ l2 


ll"h-« u ll ffi 




e 






0.75 


0.239881452 


1.034631013 


6.434091356 


0.5 


0.354575803 


1.249905596 


7.113942026 


0.1 


1.025865488 


1.738346916 


8.576177747 


0.075 


1.139434558 


1.725865966 


8.670049892 


0.05 


1.279218511 


1.67487313 


8.776573597 


0.025 


1.470238076 


1.547268561 


8.932824077 


0.0125 


1.623295808 


1.410645242 


9.070313375 


0.005 


1.793531488 


1.235799336 


9.218704868 


0.0025 


1.904725072 


1.108512517 


9.305194248 


0.00125 


2.001778889 


0.986711711 


9.371813749 


0.0005 


2.109192262 


0.835829887 


9.433204851 


0.00025 


2.176008824 


0.731445692 


9.458627896 


0.000125 


2.232164725 


0.636116593 


9.456125065 


0.00005 


2.293174219 


0.523296856 


9.360155452 



TABLE 4. Test la: Change of \\u° - ti||| w.r.t. e (h = 0.009) 



(a). u E = 20x 6 + y 6 , f e = 18000a; 4 ?/ 4 - e(7200x 2 + 360y 2 ), 

g £ = 20x 6 + y 6 , ^ = 600x 4 + 30y 4 . 

(b.) u e = a;sin(a;) + ysin(y), f £ = (2cos(x) — xsin(a;))(2cos(jy) — y * sm(y)) 

— e(xsm(x) — 4cos(a;) + ysin(y) — 4cos(?/)), 

g e = xsm(x) + ysin(y), <jf = 2cos(x) — xsin(a;) + 2cos(y) — jysin(jy). 

After recording the error, we divided each norm by the power of h expected to 
be the convergence rate by the analysis in the previous section. As seen by Tables 
[6] and [8] the error converges quicker than anticipated in all the norms. 

Test 3. In this section, we fix a relation between e and h to determine a "best" 
choice for h in terms of e such that the global error u° — u e h is the same convergence 
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h 


IK - U e h \\ L 2 


IK - 


IK ~ U%\\ H 2 


0.083333333 


4.09993E-05 


0.00268815 


0.169852878 


0.05 


1.08355E-06 


9.91661E-05 


0.011700487 


0.030656967 


3.64957E-08 


5.42931E-06 


0.000986132 


0.023836565 


7.67076E-09 


1.29176E-06 


0.000312985 


0.015988237 


4.51167E-10 


1.0898E-07 


4.04141E-05 


0.012833175 


8.8807E-11 


2.43919E-08 


1.18929E-05 



TABLE 5. Test 2a: Change of \\u e - u||| w.r.t. h (s = 0.001) 



h 






IK-<JI H 2 








0.083333333 


122.4232319 


668.897675 


3522.069268 


0.05 


69.34721174 


317.3313846 


1872.077947 


0.030656967 


43.96086573 


200.4928789 


1116.396482 


0.023836565 


41.81926563 


167.8666007 


969.5028297 


0.015988237 


27.01059961 


104.3140517 


618.4873284 


0.012833175 


19.88119861 


70.07682598 


438.4809442 



TABLE 6. Test 2a: Change of \\u e - u" h \\ w.r.t. h (e = 0.001) 



h 


IK - uIWlz 


IK ~ ul\\ H i 


\\u E ~ u c h \\ H 2 


0.083333333 


2.10771E-08 


4.91457E-07 


4.165E-05 


0.05 


5.17295E-10 


1.90347E-08 


2.72117E-06 


0.030656967 


1.77011E-11 


1.04548E-09 


2.39861E-07 


0.023836565 


3.59463E-12 


2.62405E-10 


7.50769E-08 


0.015988237 


2.03078E-13 


2.23902E-11 


9.51713E-09 


0.012833175 


3.64151E-14 


4.95624E-12 


2.69794E-09 



TABLE 7. Test 2b: Change of - u € h \\ w.r.t. h {e = 0.001) 



h 


n^-^tj L 2 


IK-i'hllffi 










0.083333333 


0.062935746 


0.122290283 


0.863654524 


0.05 


0.033106867 


0.06091104 


0.435387754 


0.030656967 


0.021321831 


0.038607272 


0.271545609 


0.023836565 


0.019597137 


0.034099981 


0.232558269 


0.015988237 


0.012157901 


0.021431653 


0.145647654 


0.012833175 


0.008152235 


0.014239078 


0.099470776 



TABLE 8. Test 2b: Change of - u € h \\ w.r.t. h {e = 0.001) 



rate as that of u° — u e . We solve problem (2.5) with the following test functions 
and parameters. 

(a) . u° = x 4 + y 2 , f = 24x 2 , g = x 4 + y 2 . 

(b) . u° = 20x 6 + y 6 , / = 18000xV, g = 20x 6 + y 6 . 

To see which relation gives the sought-after convergence rate, we compare the 
data with a function, y = (3x a , where a — 1 in the L 2 -case, a — ^ in the ii^-case 
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and a = \ in the H 2 -case. The constant, /3, is determined using a least squares 
fitting algorithm based on the data. 

Figures [5j|6] and [ljji] show that when h — es, \\u° — u%\\h 2 ~ 0(e?) and \\u° — 
w llli 2 ~ 0( £ )- We can conclude from the data that the relation, h = e 2 is the 
"best choice" for h in terms of e. It can also be seen from Figures [3]-|4] that when 
h = e, \\u°-u%\\ H mO(ei). 




005 0.1 0.15 0.2 0.25 03 '0 2 0.25 0.3 0.35 0.4 0.45 



FIGURE 1 . Test 3a. L 2 Error 
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